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We study a continuum model of overdamped self-propelled particles with an aligning interaction 
in two dimensions. By combining analytical and numerical work ,we map out the phase diagram for 
generic parameters. We find that the system self-organizes into two robust structures in different re- 
gions of parameter space: solitary waves of ordered swarms moving through a low density disordered 
background, and stationary asters. The self-regulating nature of the flow yields phase separation, 
ubiquitous in this class of systems, and controls the formation of solitary waves. Self-propulsion and 
the associated active convection play a crucial role in aster formation. A new result of our work is 
a phase diagram that displays these different regimes in a unified manner. 



I. INTRODUCTION 

Active fluids are composed of interacting self-propelled 
particles that individually consume energy and collec- 
tively generate large scale motion. Examples span many 
length scales, ranging from animal herds PQ, schools of 
fish [2], bird flocks [3j and insect swarms [4], to bacterial 
colonies [5HZ] and the cell cytoskeleton [8 j. In this pa- 
per, we consider the overdamped dynamics of a collection 
of self-propelled particles subject to local aligning inter- 
actions. The model is relevant to various experimental 
systems, including motility assays [9 , suspensions of cy- 
toskeletal filaments [TO] , self-chemotactic bacteria such as 
E-coli in convection-free geometries [11], and inanimate 
systems such as vibrated granular monolayers [12] and 
chemically driven nano-rod suspensions [T3] , 

The goal of the study is to identify universal hydrody- 
namic mechanisms for the emergence of complex struc- 
tures in systems exhibiting collective motility. To this 
end, we focus on a simple and generic macroscopic de- 
scription in terms of a conserved density field and a 
collective velocity or polarization field. The continuum 
model is expected to capture the behavior of the system 
on length scales large compared to the size of the indi- 
vidual units and time scales long compared to the micro- 
scopic interaction times and the frictional time scale set 
by the medium. Such a description was first considered in 
the pioneering work of Toner and Tu [T4HT6] who found 
that self-convection inherent to active particle flows sta- 
bilizes long range order in two dimensions. In this work 
we study this generic coarse-grained description analyti- 
cally and numerically, map out emergent inhomogeneous 
structures, and identify the mechanisms underlying their 
formation. 

We study the system as a function of two parameters, 
the self-propulsion speed wq of the active particles and a 
parameter A that incorporates the effect of inter-particle 
interactions. The theoretical model of Toner and Tu [14]- 
[T6] yields a mean field transition from a disordered state 
to an ordered "polar" state (i.e., a state with nonzero 



1.8 




\ 
\ 


po = 1 .05 

y 


1.6 
1.4 




1 A 
\ 




1.2 




\ 


y 


A 1 

0.8 




\ 


s 


0.6 








0.4 








0.2 

n 




H x 





A 



S 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



wo 

FIG. 1. (Color online) Left Panel : Phase diagram in 
the (A,u>o) plane for po = 1.05. The dashed (blue online) 
and dashed-dotted (red online) lines are the neutral stability 
curves for the L (Eq. (5)) and T (Eq. Q) modes, respec- 
tively. The shaded regions describe the stable states obtained 
numerically: a homogenous polar state (H), a regime of prop- 
agating solitary stripes (S), a regime of stationary asters (A), 
and moving localized polar clusters (B). The domain size is 
(in dimensionless units) 128 x 128 and the equations were 
integrated up to 5 x 10 4 diffusion times. The moving polar 
clusters are a result of finite system size and the finite time of 
integration. Right Panel : Snapshot of a propagating stripe 
(bottom) and a stationary aster (top), the color scheme de- 
notes density (increasing from blue to red). 



mean velocity) that is controlled by the density of the 
active units. This density is not, however, an external 
control parameter as in conventional equilibrium phase 
transitions, but rather is convected by the order param- 
eter. This coupling renders the dynamics of the system 
self regulating, in that the state of the system is deter- 
mined by the interplay between particle convection and 
tendency to local alignment, rather than by an externally 
controlled density of active particles. 

Our main result is the phase diagram shown in Fig. [T] 
that displays the various dynamical states obtained by 
varying wq and A, for a fixed density above the mean- 
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field order- disorder transition. As expected, the homo- 
geneous polar state (H) is unstable in most of the param- 
eter space. It is replaced by one of two robust structures 
in different regions of parameters. The first is a state of 
propagating solitary waves consisting of high density or- 
dered swarms moving through a low density disordered 
background (see Fig. fl Region (S)). The second is a sta- 
tionary aster (see Fig. I] Region (A)). 

Density waves of the type found here are ubiquitous 
in bacterial systems [TTJ [T7H22] and have also been ob- 
served in dense motility assays of short actin filaments [9 . 
Theoretical studies showing the emergence of travel- 
ing wave structures have included diffusion models with 
chemotactic gradients [23H26] , numerical simulations of 
agent based models, such as the Vicsek model [27] and 
coarse grained theories [28H3Q] . In particular, Bert in et 
al [29] have recently pointed out that the traveling den- 
sity stripes are solitary waves, rather than a nonequilib- 
rium pattern of the system. In this work, we show that 
these density waves are an inevitable consequence of the 
self-regulating nature of self-propelled particle flows and 
can be viewed as a coexistence between two stable phases 
of the system, namely a high density ordered swarming 
state and a low density disordered state. 

Asters are ubiquitous in cell biology in processes such 
as the formation of the mitotic spindle [3TJ [32]. They 
also occur in in-vitro suspensions of cytoskeletal filaments 
and motor proteins [10] [33] [34] . Theoretical models of 
mixtures of cytoskeletal filament and motor proteins do 
indeed yield aster formation [35H4T] , arising from the 
dynamics of the motor proteins and/or the flow of the 
solvent. Asters have not, however, been obtained before 
in Toner-Tu continuum models of self-propelled particles 
because they only occur for larger value of the effective 
nonlinearities than considered in previous work. In addi- 
tion, our work identifies a universal hydrodynamic mech- 
anism for aster formation in the change in sign of the 
effective nonequilibrium compressibility of the system, 
combined with the active self convection. 

The layout of the paper is as follows. First, we review 
the hydrodynamic theory and describe the key features 
that control the emergent structures. Then, we carry out 
a linear stability analysis of the homogeneous swarming 
state. Next, we report the results of a numerical solu- 
tion of the nonlinear deterministic equations and discuss 
the mechanisms underlying the formation of the propa- 
gating density waves and stationary asters. Finally, we 
conclude with a discussion aimed at placing this work in 
the context of the existing vast literature on active polar 
fluids. 



II. THE CONTINUUM MODEL 

We model the overdamped dynamics of a collection of 
self-propelled particles. The only conserved field is the 
number density p(r, t) of active particles. In addition, to 
describe the possibility of states with polar orientational 



order or collective motility, we consider the dynamics of 
a vector field, r(r,t) = p(r, t)P(r, £), that represents a 
polarization density. Here, P(r, t) is an order parame- 
ter for polar orientational order. Its magnitude measures 
the amount of orientational order and its direction rep- 
resents the Goldstone mode associated with the spon- 
taneously broken rotational symmetry in the swarming 
state. The dynamical equations associated with these 
quantities were first constructed phenomenologically by 
Toner and Tu [14] and have more recently been derived 
from various microscopic models [28, 29 , 38 , 40, 42 , 43 . 
They are given by 



dtp 



- V • (w r - DVp) , 



(1) 



9 t r+Ai(r • V)r = - [a 2 (p) + a 4 \r\ 2 } r + KV 2 r 

-^ 1 Vp+^V|r| 2 + A 2 r(V-r). (2) 

The density equation, Eq. [T] is a conservation law with 
a mass flux controlled by the sum of convection of the 
active particles at the self-propulsion speed wq and dif- 
fusion at rate D. The structure of the polarization equa- 
tion, Eq. [2] reflects the fact that r plays a dual role: on 
one hand r/p is the orientational order parameter of the 
system, on the other w^r/p represents the mean flow ve- 
locity. The first two terms on the right hand side of Eq. [2] 
control the mean-field continuous order-disorder transi- 
tion to a state of collective motility, with a 2 (p) a param- 
eter that changes sign at a characteristic density p c , and 
a±(p) > 0. The term proportional to Ai describes self- 
convection. It is the analog of the finite Reynolds num- 
ber convective nonlinearity in the Navier- Stokes equa- 
tion. Unlike the fluid, this overdamped system does not 
possess Galileian invariance as the particles are moving 
relative to a medium. As a result, Ai 7^ 1/p. The lack 
of Galileian invariance also allows other convective terms 
proportional to A2 and A3 appearing on the right hand 
side of the equation. In contrast to the term proportional 
to Ai, these terms also have an equilibrium interpretation 
and can be present in an equilibrium polar or ferroelec- 
tric fluid [44]. The term proportional to w± is essentially 
a pressure gradient and is unique to systems with polar 
symmetry. Before proceeding further we point out the 
two crucial properties of the above equations that con- 
trol the formation of emergent structures in this system. 

Dynamical Self-regulation. It is useful to compare 
Eq. [2] to that for the order parameter of an equilibrium 
lyotropic polar liquid crystal such as a smectic C Lang- 
muir monolayer [45]. The important difference is that in 
the equilibrium system the density controlling the order- 
disorder transition is externally tuned. In the case of a 
collection of self-propelled particles, the density is con- 
verted by the order parameter itself through the wo term 
m Eq. [1] In this sense the amount of order is itself reg- 
ulated by the dynamics of the system and this coupling 
is rendered non-local by the convective terms. This is a 
crucial feature of the dynamics of self-propelled systems 
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and plays a central role in controlling the formation of 
emergent structures. 

Negative effective compressibility. Two terms on 
the right hand side of Eq.[2]that are along the direction of 
the spatial gradient, namely w{¥ p — 4^- V|t| 2 represent 
pressure gradients, with w\p — 4f |t| 2 the effective pres- 
sure. The first term is the ideal gas part of the pressure. 
The second term arises due to the intrinsic tendency of 
polar systems to splay [44]. For A3 > 0, as obtained from 
all microscopic derivations [20 [30], it describes a lowering 
of the pressure due to ordering of the velocities. When A3 
becomes large the system develops an effective negative 
compressibility. This is the central property that controls 
the physics of this system in the interaction dominated 
regime. 

Given the large number of parameters in the hydrody- 
namic equations, Eq. [T} [2] we simplify them as follows. 
We choose a 2 = v (1 — p/p c ) an d «4 = [y / 'p 2 ) (1 + p/p c ), 
with v a characteristic kinetic frequency. This yields a 
continuous (in dynamical systems terms - a supercritical 
pitchfork) mean-field phase transition from an isotropic 
(r = 0) to a homogeneous, polar or swarming state 
(|r| > 0) at the critical density p = p c and ensures that 
\ r \/Po —> 1 for p ^> p c - We further assume D = K. 
In the following we measure time in units of v~ x and 
lengths in units of (D/v)*. Without loss of generality, 
we also set the critical density p c = 1. To reduce the 
number of independent parameters the numerical work 
is carried out for w\ — wq and Ai = A2 = A3 = A. The 
equations then involve three parameters: (1) the mean 
density of the system, po, which determines the distance 
from the order-disorder transition, (2) the convective ve- 
locity u>o, and (3) A which is controlled by the strength of 
interparticle interactions. The hydrodynamic equations 
in dimensionless simplified variables are given by 

d t p = - V • (w r - Vp) , (3a) 
d t r -(a 2 + a 4 r 2 )r - w Vp + V 2 r 

+A (r a Vr a + rV • r - r • Vr) , (3b) 

where wq and A are now dimensionless. 

III. LINEAR STABILITY ANALYSIS 

In this section, we examine the linear stability of the 
homogeneous solutions to Eqs. ([!]) and We use di- 
mensionless variables, but to highlight the role of each 
term in hydrodynamic equations, we retain the differ- 
ent parameters Wi and A^. The continuum equations 
have two homogeneous, stationary solutions, both with 
P — Po — constant: an isotropic state with r = for 
po < p c = 1 and a polar state with r 7^ for po > p c . 
We focus here on the polar (or collective motility) state. 
Without loss of generality, we take the direction of polar- 
ization to be the x axis of our coordinate system. The ho- 
mogeneous polar state is then characterized by r = tqx 



with r = Poy/(po- Pc) I (A) + Pc)- 

We now examine the stability of this state to small 
amplitude perturbations by letting p = po + £p (r, £), 
r = To + St (r, t) x + 5r± (r, t) y . Introducing the Fourier 
representation x(q, t) = J dre lCL ' r x (r, t), the linearized 
equations are given by 

d t Sp = -q 2 Sp + iwi(q\\ St + q± 5r±) , (4a) 

d t 6r = (ai + iwiq\\) Sp- (a 2 + iXr q\\ + q 2 ) Sr 

-iX 2 r q± St ± , (4b) 

dtSr^ = iw q± Sp- zA 3 t <2± St + (iXiT q\\ - q 2 ) Sr±_ , 

(4c) 

with A = A3 + A2 — Ai, q\\ — gcos#, q±_ — qsin6 
and the angle between the wavevector q and the di- 
rection of broken symmetry, x. Also, we have defined 

oti = -ro (7^ - and a 2 = -2a°, with a t > 

and OL 2 > for po > p c . We seek solutions of the form 
Sp,Sr ~ e s «(q)^ The homogeneous state will then be 
stable provided Re[s a (q)] < for all eigenvalues. The 
linear stability of the homogeneous polar state has been 
discussed elsewhere [29] [30] and will be only summa- 
rized here with focus on the aspects relevant for emer- 
gent structures. The physics is highlighted by examining 
the special cases of wavevector q along the direction of 
broken symmetry (0 = 0) and normal to it (9 = tt/2). 

Convection mediated density instability. When 
= and q = q\\, the fluctuations Sr± decouple and are 
always stable. The coupled equations for the density and 
magnitude fluctuations give the dispersion relations of 
the other two modes. In the long wavelength limit q —> 
one finds that these modes are unstable for 



OL2 



— Oi\ 
XtoWo^t < . 

ai 



(5) 



This instability is intimately related to the self regulating 
nature of the dynamics. It arises from the density de- 
pendence of the local tendency of the system to build up 
polar order, a 2 (p)^ and the fact that this density depen- 
dence is rendered nonlocal by the convective coupling of 
density to r proportional to wo and w\. The location of 
the longitudinal instability line defined by Eq. ^ in the 
(A, wq) plane is shown by the dashed (blue online) line in 
Fig. [l] for a fixed value of density. Alternatively, one can 
identify a density p^(A, wq) such that magnitude fluctua- 
tions in the order parameter destabilize the homogeneous 
polar state for pf> > p c , leading to the onset and growth 
of density inhomogeneities.This longitudinal instability 
has been discussed previously in Refs. [29 and [30] and 
will be referred to here as the convection mediated den- 
sity instability. We want to emphasize two features of 
this instability relevant for this work : i) this instability 
is mainly controlled by the third term in Eq. (5| and as 
such is model independent and arises purely due to the 
convective coupling between the collective velocity and 
density; ii) as the order-disorder transition is approached 
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from above ot^ —> and wq c ~ (po — Pc) 1 ^ 2 , i-e., near p c 
the mean-field ordered state becomes unstable for van- 
ishingly small wq. 

Splay induced negative compressibility. We now 

consider the dynamics of fluctuations in the direction or- 
thogonal to the polarization, i.e., for = | or q — q±. 
In this case the modes are all stable and diffusive near 
the isotropic-polar mean- field transition [30]. Far from 
the transition fluctuations in 5r decay rapidly and can 
be eliminated in favor of Sp and 5r±_ by letting St « 

(S) ^ " (^r) Sr± ' Substituting this in Eqs. (|4a) 



and (4c) we find that St± fluctuations, which in this case 



describe splay deformation of the order parameter, are 
unstable for 



wi - A 3 t — < . 

OL2 



(6) 



This splay instability is controlled by nonlinear couplings 
proportional to A3 and occurs when the effective com- 
pressibility of the system as determined by the third and 
fourth terms on the right hand side of Eq. ([2| changes 
sign. The parameter A3 is in turn determined by interac- 
tions in the system [29j [30] , hence the precise location of 
the instability line is model dependent. 

To summarize, we have identified two mechanisms that 
render the system linearly unstable and lead to growth of 
density fluctuations. In the next section, we go beyond 
the linear stability analysis and use detailed numerical 
computations to characterize the spatially (and possibly 
temporally) inhomogeneous states that replace the homo- 
geneous solution in the unstable region of parameters. 



IV. EMERGENT STRUCTURES 



To study the emergent structures arising from the non- 
linear model, Eqs. (3a) and (3b) were solved numeri- 



cally in two spatial dimensions using an explicit (FTCS) 
scheme and a semi-implicit Fourier- Galerkin scheme. In 
most calculations the system was started at t = in 
a homogenous, non-polar (disordered) state, with small 
amplitude, uniformly distributed, zero-mean noise. Lo- 
cal initial density perturbations were chosen to be less 
than 3% of the mean density. 

The results are summarized in the phase diagram 
shown in Fig [I] and discussed briefly in the introduction. 
For large values of activity the homogenous polar state 
(H) is unstable and two steady, inhomogeneous states are 
obtained: (i) propagating stripes comprised of ordered 
swarms moving through a disordered background, when 
active convection exceeds the strength of nonlinearities 
wq/X ^> 1, and (ii) a stationary aster when non-linear 
effects dominate convection X/wq ^> 1. In this section, 
we discuss the mechanisms underlying the formation of 
these two emergent structures. 
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FIG. 2. Profiles of p, the density (pink in online version) 
and t x and t v the polarization fields (blue and green online) 
are plotted for a striped state when wo = 0.4, A = and 
po = 1.05. The system size is 1024 x 32 and results shown are 
obtained after integrating to 10 4 diffusion times. The arrow 
on the top right (red online) denotes the direction of motion 
of the stripes The vertical lines demarcate the trailing (T) and 
leading (L) edges of the stripe. The inset shows the details of 
the profile of density (solid line, blue online) and polarization 
(dashed line, green online) for one stripe. 



A. Solitary waves 

We begin by considering the dynamics of the system in 
the convection dominated regime. For values of wo above 
the critical value for the onset of the convective density 
instability, defined by Eq. ([5|, an initially homogeneous 
state develops hot-spots of high density that are then 
convected throughout the system due to the coupling be- 
tween r and p. These high density regions reorganize 
and grow in the direction lateral to their motion due 
to diffusion, resulting in the formation of high density, 
highly ordered planar stripes propagating at a speed of 
order wq through a low density, disordered background. 
It is important to stress the properties of these stripes: i) 
they are not a pattern in that the width and spacing of 
the stripes is not fixed but rather determined by initial 
conditions and domain size - they are instead solitonic 
wavefronts; ii) the propagation speed c is of the order 
of wq and weakly dependent on the base density of the 
system; iii) The bands are bounded by two sharp fronts - 
a leading, narrow boundary layer and a wider and more 
slowly decaying trailing boundary layer (see Fig.[2|. This 
phenomenology is the same as that found in |29f . 

Mechanism. The propagating stripe state exists 
even when we set A to zero. In addition, the onset 
of the striped state follows closely the neutral stabil- 
ity line associated with the convection mediated den- 
sity instability determined by Eq. [5j Thus, a minimum 
dynamical model for the emergence of this structure is 
given by the two coupled equations dtp = —wq8 x t and 
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FIG. 3. Illustration of the phase separation in the stripe state. 
Left Panel: Computational results obtained at T = 10 4 for 
a 128 x 128 domain with A = 0. The shaded region indi- 
cates the domain in parameter space when the polar state is 
unstable due to the convection mediated density instability. 
The density in the propagating stripe ph and the low density 
background pi for various values of the base density and wo 
are shown. The lines are a guide to the eye. Right Panel: 
Plot of the polarization th as a function of wo from the nu- 
merical simulation. The solid line is the analytical prediction 
pfcOo)- 1 
p^Oo)+i " 



d t r = —(a2+a4 ] jj2)T — wod x p. The two coupled equa- 
tions are effectively equivalent to a wave equation where 
the homogeneous nonlinear terms from the polarization 
equation provide the dispersion that generates the soli- 
tary wave structure. After transforming to a co- moving 
frame, the equations can be reduced to quadrature to de- 
termine the speed of propagation and the profile of the 
wavefront. Such an analysis has already been reported 
in [29]. Here, we present a complementary point of view. 
As stated earlier, the convection mediated density insta- 
bility occurs in the vicinity of p c . The polar ordered 
state is re- stabilized at higher densities. In this unstable 
regime, the dynamics of the system essentially yields a 
phase separation into a high density ordered state and 
a low density disordered state, both of which are now 
stable (see Fig. [3|. The degree of order in the stripe is 
precisely r h = ((p h - l)/(ph + l)) 5 , i.e., the value pre- 
dicted by the mean field theory for a state of density ph- 
The propagating nature of the ordered state results in 
the robust concentration waves observed in the numeri- 
cal solution. Phase separation is also observed in active 
nematics [46 which are also self-regulating in nature, al- 
though via mechanism different from polar convection. 

Experimental Realizations. Propagating concen- 
trations waves have been observed in dense actin motil- 
ity assays [9 and in self-chemotactic bacterial suspen- 
sions [11 . For the actin system, the numerical model- 
ing described in Ref. [9 yielded the conclusion that lo- 
cal polar aligning interactions among the filaments are 
necessary for the emergence of the propagating waves. 



These interactions need not be medium mediated and 
could arise from the fact that aligned actin filaments have 
very anisotropic friction constants for sliding along the 
direction of alignment when compared to sliding against 
it [47]. This conclusion indicates that the propagating 
waves seen in actin motility assays are indeed the prop- 
agating stripes obtained in the present model of polar 
self-propelled particles due to the dynamic self-regulation 
of the collective motility. Pattern formation in bacterial 
suspensions, in contrast, are controlled by the presence of 
gradients of signaling chemicals and nutrients and have 
been traditionally modeled in terms of coupled reaction- 
diffusion equations. Recent work has shown that some 
patterns in chemotactic bacteria can be explained by 
accounting for the fact that the bacterial motility de- 
pends on density [48 . In the present context, we want 
to focus on the so-called Keller-Segal bands, extensively 
studied in the bacteria literature. They arise from self- 
chemotaxis and are the predominant emergent structure 
in convection free geometries [11 . The alignment inter- 
action among bacteria in these systems scales with the 
concentration of the chemoattractant, which is in turn 
proportional to the concentration of bacteria. In this re- 
spect the alignment interaction is similar to the model 
considered here, where the local ordering tendency is 
controlled by that depends on the concentration of 
self-propelled particles. There is, however, an important 
difference between the two systems in that in the exper- 
iment of Ref. [11] and related studies, the direction of 
propagation of the bands is set by the nutrient gradient 
in the microchannel, hence the rotational symmetry is 
externally broken. In the model considered here, in con- 
trast, the symmetry is broken spontaneously. However, 
we note that the fact that a symmetry is spontaneously 
broken is encoded in the dynamics of the associated Gold- 
stone mode and the analysis here shows that the Gold- 
stone mode plays no role in the dynamics that give rise to 
the concentration waves. Hence, the Keller- Segal bands 
are also another realization of dynamic self-regulation in 
self-propelled particle flows. 



B. Stationary Asters 

We now focus on the other inhomogeneous stationary 
state obtained in the linearly unstable region of param- 
eters, namely a stationary aster. For X/wq ^> 1 the nu- 
merical solution of the nonlinear equations yields asters 
or —1 topological defects with radially symmetric profiles 
of both density and orientational order. Other authors 
have obtained stationary asters in models with uniform 
density [35j [36], where the aster is strictly a defect in 
the order parameter. Here in contrast, it is also a re- 
gion of high concentration. A typical aster is depicted in 
Fig. [4| The density field p is a radially symmetric func- 
tion of r = |r| about the center of a —1 defect located for 
convenience at the origin of the coordinate system. The 
density is a maximum at r = and decays exponentially 
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Radial position 



FIG. 4. Radial density and |r| 2 profiles of an aster for 
po = 1.07, A = 0.8 and wq — 0.05. The system size is 128 x 
128 and the equations have been integrated up to 10 4 diffusion 
times - these profiles remained constant even after integrat- 
ing to 5 x 10 4 diffusion times. The density field p(r) exhibits 
a point of inflection; at the same radial position the order 
parameter \t\ attains a maximum value r m ax- The aster is 
characterized by a core of size £ co and a region where both 
density and polarization decay over a characteristic length 
scale £oo. The dashed curves (purple online) are the theo- 
retical result obtained by the asymptotic analysis based on 
Eq. (Trl) and described in the text. 



far from the core to a value slightly below the critical 
density for the onset of ordering, p c = 1. Unlike the den- 
sity field, the order parameter, r(r) is non-monotonic. 
Starting with a zero value at the core, the order parame- 
ter increases almost linearly to a intermediate maximum 
value r max and decays exponentially to zero far from the 
aster's core. The point of maximum r also corresponds 
to a point of inflection in the density profile. 

We characterize the aster by two length scales, the size 
£ co of the aster core as defined as the distance from the 
aster's center to the point where |r| reaches its maxi- 
mum and the density has an inflection point, and the 
length scale characterizing the exponential decay of 
both density and polarization. This second length can be 
thought of as the size of the aster. Both length scales are 
only weakly dependent on the interaction strength A for 
fixed wq. For fixed A we find that both £ co and £00 de- 
crease as a function of wq while the maximum density and 
the maximum r increase. In other words, for fixed inter- 
action strength larger self-propulsion speed yields denser, 
tighter asters. The behavior of density and order pa- 
rameter far form the center of the aster can be obtained 
from a simple asymptotic analysis of the dynamical equa- 
tions as the nonlinearities do not play a significant role 
in determining the asymptotic behavior. Far from the 
core, assuming a radial dependence of both the density 
and the order parameter, and using the fact that r — >> 
at infinity, we can estimate the magnitude of the polar- 



ization by considering the steady state equation to lin- 
ear order in small deviations from the asymptotic values 
(p,r) = (poo<l,0) 

r 2 g+r^-[r 2 (l- Poo+ ^) + l]r = 0. (7) 

The solution to the above equation is a Bessel's function 
of the second kind, that decays exponentially as e £ ~ 
with the length scale ~ (1 — p^ + Wq)~z , consistent 
with the trends identified in the numerical work. The 
behavior near the core is determined by a complicated 
interplay between the nonlinearities and cannot be in- 
vestigated analytically. 

Mechanism. The basic mechanism underlying aster 
formation is the splay-induced negative compressibility 
discussed in the context of the linear stability analysis. 
We stress, however, that the self-regulating nature of the 
flow is critical for this instability as well in that if we set 
Wq to zero, we do not obtain any emergent structures. 
The system forms asters because this is the only struc- 
ture that can accommodate both the tendency to splay 
and that to phase separate. Unlike the solitary wave 
discussed in the previous section, the aster state is not 
universal and depends on the values of the parameters in 
Eq. [2] The self-convection term (Ai in Eq. (|2|) is criti- 
cal for the formation of this structure. When we set this 
term to zero we find that the system develops streamers 
instead, a phase separated polar ordered state where, in 
contrast to the solitary wave state, the polarization is 
orthogonal to the direction of the spatial gradients (see 
Fig. |5|. 

Experimental Realizations. Aster formation has 
been seen in purified cell extracts of microtubules and 
associated motor proteins, such as those studied in 
Ref. [lOj [33j [34]. In these controlled in vitro systems, 
with known concentrations of only a few types of motor 
proteins, aster formation can be understood theoretically 
using both continuum models [35 , 36 , 38 and simulations 
[4T] . Understanding the origin of aster- type structures in 
vivo, such as the formation and splitting of the mitotic 
spindle upon cell division, is much more challenging as in 
this case the process is controlled by a variety of compet- 
ing nonequilibrium mechanisms, including microtubule 
polymerization and the on/off dynamics of many differ- 
ent motor types [31] . Our continuum theory where aster 
formation is controlled by only two competing parame- 
ters A and wo may then provide a useful guidance for 
the modeling and interpretation of these complex exper- 
iments. In particular, one can imagine fitting the spindle 
structure obtained in experiments on cells where various 
proteins have been systematically suppressed [3TJ [32] to 
our continuum model to identify which proteins directly 
affect the model's parameters and thereby back out the 
role played by the protein in the formation of the struc- 
ture. 
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FIG. 5. Steady inhomogeneous phase- separated structures, 
we term streamers, obtained at fixed domain size 128 x 128 
when Ai = 0, A 2 = A 3 = 0.8, w 0.1 and p 1.05. The 
high density (red) region is also highly polar and ordered with 
the direction of polarity along the neutral direction. In (b) 
we quantify this by plotting the two components of r as a 
function of y. Note that r x (green dashed line) is a single 
valued function of y while r y (blue solid line) changes sign as 
we traverse the phase separated ordered region from one side 
to the other. 



V. DISCUSSION 

In this work, we have considered a generic continuum 
theory of self-propelled particles in the overdamped limit 
in two dimensions. An important new point of our work 
is the identification of dynamical self-regulation - the fact 
that the density that controls the amount of order is itself 
convected by the order parameter - as a crucial mecha- 
nism for the formation of emergent structures in these 
systems. We identify two robust structures in different 
regions of parameters: traveling solitary waves and sta- 
tionary asters. We characterize these structures using 



numerical solutions of the deterministic nonlinear equa- 
tions and identify the underlying hydrodynamic mecha- 
nisms associated with their emergence. 

The primary limitation of this study is that our sys- 
tem is overdamped or "dry" . Momentum is not conserved 
and there is no coupling between polarization and actual 
flow. The "dry" system should be contrasted with an 
active polar suspension, consisting of active particles in 
a fluid that mediates hydrodynamic interactions. In the 
suspension the total momentum is conserved and impor- 
tant nonequilibrium effects arise from the coupling of the 
associated flow velocity and local orientational order, as 
shown in recent work [49, 50 . A second limitation is the 
one elastic constant approximation that identifies the en- 
ergy cost for bend and splay. Theoretical work on models 
of suspensions of cytoskeletal filaments has shown the im- 
portance of retaining different elastic constants for bend 
and splay to understand the emergence of vortices and 
spirals [37j [ST] . This physics is missing in the present 
model. 

Finally, we note that aster formation has also been ob- 
served in related continuum models when the sign of pres- 
sure gradients tends to favor the formation of high den- 
sity regions [38, 39 . In contrast, asters were not found in 
Refs. [29j|30], even though the continuum equations con- 
sidered there have the same structure as those analyzed 
here. The reason for this is that Refs. [29j [30] study the 
continuum model obtained by systematic coarse-graining 
of specific microscopic models (the Vicsek model and a 
collection of self-propelled hard rods, respectively) and 
use the microscopic expressions for the parameters ob- 
tained in such derivations. In both cases the microscopic 
calculation yields A ~ Wq. In other words A is not an 
independent parameter and cannot be tuned to the large 
values required to obtain asters. 

The simplicity and generic nature of the theory con- 
sidered in this study has enabled us to highlight the role 
of dynamic self regulation and to show that the mecha- 
nism is universal and does not depend on the microscopic 
physics, in contrast to closely related albeit system spe- 
cific studies such as those in Refs. [39) [52]. 
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